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_ Abstract 



Radio astronomical observations have very poor signal to noise ratios, unlike in other disciplines. On the other hand, it is 
possible to observe the object of interest for long time intervals as well as using a wider bandwidth. Traditionally, by averaging 
O |. in time and in frequency, it has been possible to improve the signal to noise ratio of astronomical observations to improve the 

D ■ dynamic range. This is possible due to the inherent assumption that the object of interest in the sky is invariant over time and the 

I frequency range of observation. However, in reality this assumption does not hold, due to intrinsic variation of the sky as well as 

due to errors generated by the instrument. In this paper, we shall discuss an alternative to averaging of images, without ignoring 
subtle changes in the observed data over time and frequency, using subspace decomposition. By separation of data to signal and 
noise subspaces, not only would this improve the quality of the data, but also enable us to detect faint artifacts due to calibration 
' I errors, interference etc. 



I. Introduction 

Astronomical observations have very weak signal to noise ratios (SNR). In radio astronomy (interferometry), this weak SNR 



O 
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^ ■ is improved mainly by prolonging the observation time. This not only improves the SNR, but also sampling of the uv plane. 
A typical observation has a bandwidth of a few MHz, divided into smaller narrowband channels. Normally the bandwidth 
of a single channel is about a few kHz. After interference mitigation, calibration and imaging (also deconvolution) each of 
these narrow band channels, it is normal procedure to average them, on the image plane. The resultant averaged image has 
^ an improved SNR compared to images made by each of the narrowband channels, assuming all systematic errors have been 
OO removed. This is the typical process by which current radio telescopes like the Westerbork Synthesis Radio Telescope (WSRT), 
is able to achieve dynamic ranges in the order of 1 ,000,000 to 1 . 

However, there is a significant drawback in this averaging procedure, because it is based on a erroneous assumption: that 
the source as well as the effects of the instrument remains invariant in all averaged data. Only the noise varies and averaging 
significantly decreases the noise power. The invariance assumption is accurate to some extent, but closer scrutiny reveals that 
there is significant amount of variation: 
. • There is intrinsic variation of the sky (or the sources) with frequency (non zero spectral indices). 
The uv points scale with frequency, thus the point spread function (PSF) change in the images. 
Transient sources as well as radio frequency interference (RFI) create temporal variation. 
Calibration errors, small as they might be, create subtle variations in time as well as in frequency. 
^ ■ Indeed, by averaging, both noise and the aforementioned effects get suppressed. However, it is to our advantage that we 
detect such variations in our images. For instance, we could detect trace amounts of RFI and calibration errors etc., that appear 
only in certain time or frequency ranges. This could also enable us to discard certain parts of the data, where such effects are 
dominant, and thus improve the quality of our averaged, final image. 

Subspace techniques have previously being used in image compression [1], [2], as well as image comparison. However, to 
the best of our knowledge, this has not been used in (radio) astronomical image processing. 

The rest of the paper is organized as follows: In the next section, we give a formal representation of images and subspace 
decomposition. Next, we investigate various uses of the image subspace decomposition. Later, we consider possible extensions 
of the proposed method, by taking into account the data weights as well as the affect of PSF. Finally, we give our conclusions. 
Notation: All bold lowercase letters represent a column vector, all bold uppercase letters represent a matrix. 

^This work was first presented in the SKA calibration and imaging workshop, Cape Town, SA, December 2006 and in URSI General Assembly, Chicago, 
US, August 2008 



II. Mathematical Preliminaries 

Let us represent an arbitrary image, i.e., a set of pixels with arbitrary arrangement in two dimensions, as a vector b. Let the 
number of chosen pixels be N, so the dimension of vector b is x 1. Instead of one such image, we have a set of images, 
of the same part of the sky. Let us denote this set by X. Let us assume we have M such images, which could differ in terms 
of any of the following: 

• Observing frequency, one image per each channel. 

• Observing time, different observation dates. 

• Observing instrument, using different radio telescopes, i.e. WSRT or the Very Large Array (VLA). 

• Integration time, using snapshots 

In other words, the images in the set X have one or more forms of diversity. Let hj represent the j-th such image from X. 
We assume all the images in X could be represented by an arbitrary orthonormal basis V (N x N matrix), as in ([T]). Note that 
V is used to represent the components visible in the image, i.e., the signal, not the noise. 

b, =Va,+n, je[l...M] (1) 

In ([T]), a.j represent the components of the basis V used to represent the j-th image. The noise, which is assumed to be 
independent in each image in X is represented by n (A^ x 1 vector). Moreover, the signal or the artifacts are represented by 
Vsij in O. 

We make the following assumptions: 

• The variations in all the images in X is small enough to be represented by at most A' orthonormal vectors. 

• The noise in each image a stationary random process with autocorrelation cr^I, and is uncorrected between pixels as well 
as with the signal. 

• The noise power is much lower compared to the signal or artifacts. 
Let us consider the autocorrelation R of the images, given by (|2]). 

R = E{hh^} (2) 
Using M images, the estimate of the autocorrelation R is given by (O. 

1 ^ 

We exploit the assumption that noise is uncorrected with signal, i.e., £^{na^} = E{a.n^} = 0. Because VV^ = I, and 
we get 

R = V£;{aa^}V^ + a^I = V (^^{aa^} + a^l) V^. (4) 
From Q, we see that the column space of R is a subset of the column space V. 

co/(R) C col{\) (5) 
Using the Singular Value Decomposition (SVD) we can find an equivalent column space R = UXIU^ as in ([6|. 

co/(U) C col(V) (6) 

The singular values in 5] will characterize the signal and noise subspace eigenmodes. Because the noise has much lower 
power compared with the signal, the eigenmodes corresponding to the dominant singular values in will represent the subspace 
of V. 

Hence, we can divide eigenvectors as noise or signal according to their singular values (|7]), where Ug and Un represent the 
signal and noise subspaces, respectively. 



U = 



(7) 



Given an image bj, we can separate the signal hsj and noise hnj as in ([8]). 

b,, =U,Ufb, b,, = (I - U,Uf )b, (8) 

Note in practice, in ([3]), we normally have fewer images than the image size in pixels, i.e., N M. Hence, we would 
only have M non zero singular values in X). However, by our assumption (small enough variation), it should be possible to 
represent the signal space using fewer eigenmodes than M. 



III. Applications 

The aforementioned technique has been widely used in image compression and comparison in other discipHnes. However, 
we could find novel uses of this technique in radio astronomy. 

• Enhancement of signal to noise ratio: For each image hj in the set X, we get the enhanced version hgj in ([8]). These 
images could be used for further processing, like estimation of spectral indices of fainter sources. 

• Automatic image classification: A typical observation (say in LOFAR) consists of a few thousand channels, with an image 
being made for each one. It is assumed the majority of these images are good quality, with proper calibration and imaging. 
However, a few of these images would still have faint RFI or calibration errors. In order to automatically determine image 
quality, we could use the power of the noise component, i.e., the norm of hnj in ([5]). For a poor quality image, this should 
be much higher compared to a good quality image. 

• Continuum subtraction, spectral lines: Let us assume that the j-th image contains a weak spectral line, hidden under a 
strong continuum signal that is common to all the images in the set X. Just by looking at this image hj , it might not be 
possible to observe this spectral line because of the continuum signal. However, by subtracting this, in b^j, we should 
be able to observe this. 

• Detection of faint, narrowband RFI: The same technique applied to spectral lines could be applied to detect faint, 
narrowband RFI. It is assumed that the RFI affects only a few images in the set X. 

An implementation of this method which is capable of decomposing FITS files can be downloaded from [3]. 

IV. Example 

In this example, we have considered a WSRT observation with faint narrowband RFI, which has not been detected by the 
flagging algorithm. We have about 32 channels of data, each with a bandwidth of 100 kHz centered around 140 MHz. We 
have applied the proposed technique on this dataset. The results are given on Fig. [T] 

In order to reduce the computational cost (for an image with N pixels, the cost of computation of the SVD is 0{N^)), 
we have used a divide and conquer approach. Given a set of images with a large number of pixels, we first divide them into 
smaller sub-images, using an arbitrary boundary. For each set of these sub images (one taken from each large image), we 
apply the proposed technique. This saves the computational cost significantly. However, this could also lead to new artifacts. 
For instance, the tile like structure appearing on Fig. [T] (d) is due to these sub-images. But as seen from this image, these 
artifacts are less significant to affect the outcome. 

V. Conclusions 

We have proposed a technique to improve the quality of the data in astronomical observations using subspace decomposition. 
Application of this technique to real data has shown us that it is superior to the normally used averaging techniques. 
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Fig. 1. Partially deconvolved residual image around Cyg A using WSRT data: {(a)} residual image made by a single channel of data, {(b)} mean image 
using 32 channels, {(c)} difference image between (a) and (b), {(d)} residual remaining after removal of signal subspace from the image {(a)}. It is clearly 
evident that the image (a) contains faint RFI, seen on (d) as diagonal lines in black. These appear also on the difference image (c) but it is hard to distinguish 
them because of the artifacts remaining. The artifacts on (c) are outlined using the red polygons. By using the proposed technique, we are able to remove 
most of these artifacts. The tiled structure on (d) is due to the divide-and-conquer approach in applying this technique. 



